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1. UPDATED OBJECTIVES 
1.1 The new challenges 

Aero and space borne remote sensing have made it possible to monitor the mass and 
energy transport at various scales within the cryosphere-hydrosphere-atmosphere system. 
The recent surface mass balance (the rate of net gain of snow and ice at a geographic 
point) map for the Antarctic ice sheet is constructed by interpolating sparse in situ 
observations (about 1,800 points) with empirically calibrated satellite data of passive 
back emission of microwaves (Vaughan et al, 1999). The digital elevation model 
obtained from satellite radar altimetry is used to improve the delineation of the ice flow 
drainage basins (Bamber, 1994; Bamber and Huybrechts, 1996; Vaughan et al, 1999). As 
important as these results are, the uncertainty remains up to about ±2mm/yr. of eustatic 
sea level change with the net imbalance (Jacobs et al, 1992; Vaughan et al, 1999). In 
other words, we are still unable to determine even the sign of the contribution of the 
Antarctic ice sheet to contemporary sea level change. The problem is more likely with the 
discharge rather than accumulation. 

Transition from snowfalls to an icy layer can be regarded as a process of gravitational 
stabilization, a process may take centuries (e.g. Liboutry, 1998). Therefore, the in situ 
observations in the accumulation zones, i.e. most of the surface cover of Antarctic, can be 
progressively updated and back checked to further constrain the accuracy of mass balance 
maps (Bull, 1971; Giovinetto and Bentley, 1985; Vaughan et al, 1999). The opposite is 
true for the processes of attrition and ablation. The surge and collapse of ice streams are 
nonlinear functions of a number of dynamic variables and empirical parameters, few of 
which are well constrained. Satellite radar altimetry with the footprint of tens of 
kilometers is not adequate for monitoring the surge and collapse of ice streams. The 
technique of interferometry could recover the discharge of streams from the elevations 
observed with synthetic aperture radars (Rignot et al, 1996; Rignot, 1997; 1998), yet the 
irregular topographies, including the calving ice walls and melting ice lobes at the ice 
sheet margin, often corrupt the fringes. The highest resolution is achieved by laser 
altimetry with footprints of only tens of meters. Even with perfect altimetry, uncertainty 
in the inference of mass balances is still not well constrained without addressing the 
issues of unstable snow masks, unknown local density profiles, unknown internal flows, 
crustal responses etc. 

Basal lubrication and friction (often due to the roughness of the ice-bed interfaces) 
have strong control over the sliding and viscous deformation of the lower part of the ice 
stream. Thus, on the verge of break up, when shear stresses are released, the change in 
surface elevation may not synchronize with the flow in the lower part of the ice sheets. 
Crustal deformation in response to the ice mass redistribution, known as glacial isostasy, 
is not only an error source for the inference of mass balance with altimetry, it feeds back 
to the control mechanism for the advance or retreat of ice sheets (Hughes, 1987; Denton 
et al, 1992; Hughes, 1998). As is well known, the residual crust uplift since the 
termination of the last Quaternary cycle is a source of error for using altimetry to infer the 
change of ice sheet thickness in the Greenland and Antarctic areas. At the same time, it is 
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a major source of data for the study of mantle rheology (e.g.Peltier & Andrews 1976; 
Mitrovica & Peltier, 1993, Hager, 1989; Fang & Hager, 1996, 1999). 

A cross constraint for the ice mass imbalances is the sea level change. The 1-2 mm/yr 
sea level increase derived from tidal gauge records (Douglas, 1991, 1997; Woodworth et 
al, 1992) has been of great interest in recent years for its close relation to the publicized 
issue of global warming. But the physical significance of these figures is limited by the 
terribly sparse and uneven distribution of tidal gauge records. Even a static ocean cannot 
be perturbed uniformly (Farrell & Clark, 1976; Dahlen, 1976; Conrad & Hager, 1997; 
Fang & Hager, 1999; Mitrovica et al. 2001). Most of the published results are based on 
tidal observations along the North Atlantic coasts. It is unclear if these measured sea level 
rises apply to other coastal areas all over the world, and more uncertain of how it 
represents the entire sea level change associated with the ice mass imbalance, because 
static changes in the magnitude of l-2mm/yr can be easily over shadowed by low 
frequency ocean circulation modes (Wunsch, personal communication). 

1.2 Adjusted research strategy 

In response to the new challenges, we adjusted our main goal of this project to 
establishing new constraints for the mass balances, and more importantly the total 
imbalances, in the Greenland and Antarctic areas by quantifying their signatures from 
gravity observations of different scales, local regional, and global. 

For local analysis , it was first suggested by our group to utilize the observed crustal 
deformation to “weigh” the time variation of ice sheets (Hager, 1991). This idea was 
modified and tested by Wahr’s group (Wahr et al, 1995; Wahr & Han, 1998; Larson & 
van Dam, 2000). The basic technique in Wahr’s method is to combine the GPS vertical 
and absolute gravity measured at the same sites near the margin of ice sheets to remove 
the ancient signals from postglacial rebound to avoid it being mistaken as due to modem 
glacial change. There are two issues left by their analysis: (i) independence of the 
technique on the largely unconstrained mantle viscosity structures, (ii) the formulation for 
quantifying the ice signature in the combined observables. We solved the first problem 
during this project and published the results (Fang & Hager, 2001). Details of the results 
and current progress will be outlined in section 2. 

For regional analysis , we developed and implemented the technique of spherical 
wavelets. A wavelet transform is to use a small wave of zero mean, scaled and 
translated, to scan through the space of data. 

A unique quality of wavelets is its ability to isolate and enhance small, however small, 
and discontinuous features, thus, it is ideal for localized data sets with irregular 
boundaries and sharp contrasts such as the distribution of ice sheets, plume like structures 
in seismic tomography (e.g. Bethoux et al, 1998; Bergeron et al. 1999; Bergeron et al, 
2000). Another noticeable quality of wavelets is that it provides spatial and harmonic 
(frequency) localization simultaneously (subject to Heisenberg’s uncertainty principle). 
This is especially important for gravity analysis, which is a highly lumped signal. 

We worked on spherical wavelets in 1994 in collaboration with Dr. Mark Simons, then at 
MIT. We independently constructed what is equivalent to the locally supported spherical 
wavelet of Freenen & Windhuser (1996) employed in tracing the source of the 
conspicuous gravity low in the Hudson Bay area (Simons & Hager, 1997). The general 
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theory for spherical wavelets has been developed by Freenen & Windhuser (1996, 1997). 
Since we all use small waves (non-uniform windows) as a set of spherical caps with 
wiggled cross sections, we shall call them isotropic spherical wavelets (ISW). We 
applied wavelets to characterize the Last Quaternary ice sheets, and to quantify the 
exponential-like relaxation curves near the centers of the ice sheets. Results of these 
works are published (Fang & Hager, 2002). Our current investigation is on quantifying 
the mass balance signatures of the Greenland and Antarctic Ice sheets. The published 
results and ongoing progress will be outlined in section 3. 

For global analysis , we look at the lower degree harmonic coefficients (e.g. Cheng et 
al, 1997). These low degree coefficients are powerful constraints for the total mass 
redistribution. Yet they are highly lumped signals, even the core may have a contribution 
(Fang et al, 1996). We have found during this NASA investigation that the time varying 
harmonic coefficients, as low as degree3, are sensitive to the ice history. There have 
been a number of follow up confirmations, since we first reported the result to an IAG 
assembly in 2000. Currently, we are collaborating with the group of the Center of Space 
research at Texas University to further quantify the sensitivity. This part of the work is 
reported in section 4. 

For sea level analysis . In conjunction with previous NASA project, we have 
succeeded, during this investigation, in formulating and implementing the equilibrium sea 
level by means of least potential energy principle (Fang & Hager, 1999). The theory 
could be extended to a non-equilibrium ocean by including the non-gravitational forces 
and local tectonics to address the issue of mass budget in the ocean. But the latter idea is 
a little too ambitious. Collaborative efforts are being made by our group and the geodetic 
group in the Hong Kong Polytech. University on the analysis of tidal gauge records; 
details will be reported in section 5. 

2. PROGRESS OF LOCAL ANALYSIS 


2.1 Published results 

An increasing number of instruments, including superconducting gravimeters, absolute 
gravimeters, and GPS receivers have been operating in the Greenland and Antarctic 
regions. Crustal deformation in these areas is strongly influenced by both postglacial 
rebound and contemporary ice mass redistribution. To quantify the ice mass signature, 
we first have to get rid of the ancient signals. We seek to determine whether, when the 
viscosity profile in the Earth's interior is unknown, the effects of viscous relaxation in the 
memory of surface mass change can be separated from the effects of present day mass 
variation by combined measurements of vertical displacement and absolute gravity. The 
key issue is to establish a relationship between the vertical load Love number h n and the 
gravity load Love number k n at each harmonic degree n. We find through dynamic 
analysis (see Fang & Hager, 2001a for detail) that for an incompressible Earth, 
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where p E is the average density of the whole Earth, and p n is what we call the pseudo- 
surface density as a function of harmonic degree n. 
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The radial integration is from the center of the Earth to the surface a. 

Equation (1) is independent of rheology parameters in the interior. For PREM model 
(Dziwonski & Anderson, 1981) the average pseudo-surface density is <p>~ 3.5 . From 
(1) and (2) we come up with the analytical expression for the constant A deduced 
empirically by Wahr et al (1995): 
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Let us combine the observed vertical deformation u(t) with the incremental gravity 8g(t) 
by forming the new observable 

A{t) = u(t)~ A8g(t) (4) 


It can be shown that A(t) is primarily due to the variation of ice mass (Wahr et al. 1995; 
Fang & Hager, 2001a). Relation (4) is one of the major steps in local analysis. 

For a compressible Earth, equation (10 should be modified as (Fang & Hager, 2001a) 
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represent the elastic and viscous compressions respectively. We assert that the term of 
viscous compressibility x p is negligible for a Maxwell Earth, since it is viscously 
incompressible. This assertion is confirmed numerically as shown in Fig. 1, 2 where we 
generated three markedly different viscosity models (Fig. 1), and then calculate the 
pseudo-surface densities for each one of the models (Fig. 2). Our assertion is confirmed 


Hager & Fang 


5 


Final report 




Fig. 1 Test viscosity models. Model A 
(solid line) is a fairly uniform profiled 
advocated by Tushingham & Peltier (1991). 
We adopt it here as a reference for the two 
extremes B and C. The lithosphere for 
model C (short dashed line) is the same as 
for model A. Model C and B are generated 
using the empirical parameterization 
introduced in Fang and Hager (1999). 


Fig. 2 Pseudo-surface densities for 
incompressible and compressible Earth 
models. The compressible and 
incompressible elastic p n are displayed 
with gray lines. The thin dot line indicates 
the true surface density of the 1066B Earth 
model. The dark solid, long dashed and 
short dashed lines are the viscous 
compressible p n calculated by means of a 
least square procedure with respect to the 
time variable. 


by the convergence of the viscous pseudo-surface densities of those viscosity models (see 
details in Fang & Hager, 2001a). We plot in Fig. 3 the predicted present day change of 
gravity due to postglacial rebound in the Antarctic area for each viscosity model in Fig. 1 
along with their residuals. The present rate of residual is defined as dSg / dt = A~'du/dt 
We can see clearly from Fig. 3 that the residuals are about 5% and less in areas where the 
ice signals are strong. A better reduction of the residuals can be achieved if we replace 
the average pseudo-surface density with the true pseudo-surface density for each single 
harmonic. But this procedure is practically meaningless, since it is very difficult to obtain 
the global harmonic coefficients for the present day vertical deformation field for the 
crust. 


2.2 Results to be published 

There three issues remain to be solved before we can quantify the ice mass signature 
from the hybrid observable A(t). First, the residual control. The pseudo-surface densities 
are function of harmonic degrees as we can see from Fig. 2. It is readily shown by 
examining single harmonic degrees that for the long waveband (degree 2-9), the 
dominate waveband (degree 10-35), and the short waveband (degree >30), the average 
values for the constant A are 5.7, 6.4, 7.4 respectively. By taking the procedure (4) and 
using the constant given by (3), we could characterize the residuals of viscous effect by 
three residual values of A: -0.7, 0.1, and 0.1, respectively representing the three 
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Fig. 3 Predicted present day 
rate of change of the gravity 
field due to postglacial 
rebound in the Antarctic area 
and the residuals. Calculation 
are based on the admissible 
viscosity models displayed in 
Fig. 1, and the ice sheet model 
ICE-1A which is a modified 
version of the ice model of 
Peltier & Andrews (1976). 
Major difference between ICE- 
1 and ICE-1A happens to 
occur in the Antarctic area, 
since ICE-1 is reconstructed 
on for the northern 
hemisphere. Antarctic ice is 
first added by Nakada & 
Lambeck (1987) based on the 
Quaternary ice sheet map of 
Denten & Hughes (1986). 


waveband. It is not clear yet whether the simple average (3) is the optimal value for 
keeping the residuals of the ancient signals minimum. We plan to look at the issue by 
weighing contributions from the three wavebands. Viscosity is no longer the issue in this 
task, but ice model is. We are conducting numerical tests on all existing last Quaternary 
ice sheets models (ICE-1A, ICE-3G, ICE-4G etc) to learn the sensitivity of the optimal 
solution for A on the ice sheet model. We also plan to scale the Antarctic portion of these 
ice models artificially to see under what condition the residuals will be out of control. 
These experiments on Holocene ice models are complimentary to the viscosity control 
analysis outlined above. We expect to come up with a optimal value of the constant A 
with an upper bound for its applicability. 

The second issue concerns the compression term. Equations (1) to (6) not only 
represent the physics for removing the postglacial rebound signature, but also the 
foundation for deconvolutions for modem ice mass signatures (see Fang & Hager, 2001 
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for details). Effects from elastic compressibility have been perceived in geodetic 
investigations for Quaternary ice sheets but without formulation (Ekman, 1992; Ekman, 
& Makinen, 1996). The compression term x e n expressed by (6) can be regarded as the 
compression load Love numbers. They convolute with the ice height variation £, to give 
rise to the compressibility correction for gravity disturbance 

5g x (d,(p,t) = 2> + lKP n (cosa)^' (7) 

JV1 E n 

where P n is the Legendre polynomial, a denotes the angle between the fixed point (0, (p) 
and the moving point (0', cp'), dQ. stands for the area element on the surface of a unit 
sphere. We are investigating these compression Love numbers on general terms. Note, 
the familiar gravity kernel (1 + k n ) has to be replaced by the new kernel (1 - x e n ) after the 

removal of postglacial effects. We are conducting systematic forward modeling for the 
gravity field with both kernels. As an input of the ice mass model, we will adopt the 
existing mass balance maps for Greenland and Antarctic areas (e. g. Reeh et al. 1999; 
Zwally & Giovinetto (2000); Vaughan et al. 1999). Special attention is paid to the short 
waveband. The short waveband is not affected significantly by postglacial rebound, but it 
may have significant effects from modem ice mass balance, especially in the areas of ice 
sheet margins where most of the polar stations are located. 

The third issue that we are dealing with is about the distribution of sites. Application 
of our formulation to analyze the data collected from the polar stations is straightforward. 
Three or four sites can only serve as benchmarks for constraining the regional mapping of 
mass balances, although they are extremely important benchmarks. Taking advantage of 
the knowledge learned from the work outlined above, we are investigating the possible 
distribution with minimum number of gravity-GPS stations for a robust constraint of the 
ice mass redistribution. As far as we know, several projects are already underway for 
similar goals. What makes our method unique is that with solid physics and explicit 
formulations (1) to (6) (also see Fang & Hager, 2001 for more) we could go beyond pure 
empirical experiments in dealing with the separation of ancient secular ice signatures and 
modem short period ice signatures, and in quantifying the ice mass balances in the 
combined observables. An important application of the technique of combination is to 
combine the laser altimetry with the time varying gravity like the GRACE data (this is 
already beyond local analysis) for the estimate of ice mass balance. The basic idea of 
such combination was suggested by Wahr & Han, (1998) and tested empirically by 
Velicogna & Wahr (2000). It is relatively easy for GRACE to separate the secular term 
(presumably the postglacial rebound signal) and short period terms (here we ignore the 
uncertainties associated with the unknown continental hydrological system and 
inadequate lifespan just for the sake of discussion). It is not easy however to isolate the 
secular terms from altimetry observations, since altimetry only measures the ice surface 
rather than the basal surface. We propose to calibrate the secular term isolated from the 
GRACE gravity rates to derive the postglacial rebound correction for altimetry, by using 
the basic formulation (1), (2), (5), and (6). Here we skip Wahr’s relations (4), because 
GRACE and altimetry data could come in harmonic degrees, therefore, we could perform 
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the combinations degree by degree using the exact relation (1)! This has not been tested 
before. We believe that it will yield exciting results in terms of removing the post glacial 
rebound effect from the laser altimetry, thus, give rise to a robust constraint on the mass 
balances in Greenland and Antarctic areas. The combination technique could even 
extended to spherical wavelet analysis which we will elaborate in the next section. 

3. PROGRESS OF REGIONAL ANALYSIS 

3.1 Published results 

Ice mass in Greenland and Antarctic areas are confined to their continent beds by 
surrounding oceans. This is the main reason, aside from its unique qualities, we employ 
spherical wavelets for regional analysis. Here we demonstrate how wavelets is used to 
characterize the Laurentide and Fannoscandian Holocene ice sheets and to constrain the 
mantle viscosity structure. These are the long wavelength problems as oppose to the short 
wavelength contemporary ice mass balance problems. But the basic technique is shared 
by both problems. Since the technique is not popular yet, we include some background in 
the problem description. 

A unique feature shared by the two major Holocene ice sheets, Laurentide and 
Fennoscandia is the apparent exponential RSL curves at their central regions. This 
unique feature laid the foundation for Haskell's (1935) classic analysis to derive the 
relaxation time x for a single exponential mode exp(-£/z), at the characteristic wavelength 
of the Fennoscandia ice load. In order to achieve an exact exponential mode in a half 
space, the real ice dome has to be distorted into an exact harmonic undulation extending 
to infinity. An opposite approach to recover the apparent exponential RSL is to preserve 
the physical integrity of the ice dome by summing up all harmonics. The later approach 
has become standard in modern studies (e.g., Farrell & Clark, 1976; Peltier, 1976; 
Lambeck et al. 1990; Mitrovica, 1996; Fang & Hager, 1999). The so called decay time T 
in Mitrovica's (1996) analysis is identical with Haskell's relaxation time in the sense that 
they are obtained from fitting the same curves by essentially identical expressions 
(different only in spatial and temporal references). But there is a fundamental difference: 
a decay time is obtained by fitting a RSL curve calculated by summing up all harmonic 
degrees. 

The question here is whether the distribution and the history of retreat of an ice sheet 
can make the ice-ocean-earth system to achieve a state where the exponential like RSL at 
its center is related to a single harmonic relaxation without sacrificing the physical 
integrity of the ice dome. This is analogous to the problem Heisenberg encountered, 
which leads to his famous uncertainty principle. According to the uncertainty principle, 
the answer is no, because Haskell’s problem is to seek an exact harmonic localization and 
the modem approach (e.g., Mitrovica, 1996) is to seek an exact spatial localization. It is 
impossible to achieve both simultaneously. On the other hand, also according to the 
uncertainty principle, a compromise state could be reached at the expense of ambiguities 
in both spatial and harmonic localization. The robustness of a spatial or a harmonic 
localization depends on the physics of the system. We use ISW to quantify such a state, 
and dynamic analysis to test the robustness of harmonic localization. 


Hager & Fang 


9 



Final report 


The most appealing feature of a locally supported ISW for our purposes is its close 
relation to a smooth window. The bandpass filter is nearly symmetric about the central 
frequency (see Fig. 5 below). Mathematically, a spatial localization by an window in the 
spatial domain corresponds to a harmonic localization in the form of a band-pass filter in 
the harmonic domain. The size of the window and the width of the filter are related by 
the famous “uncertainty principle” first quantified by Heisenberg in his dealings with the 
wave-particle duality of electrons. 

To demonstrate the effectiveness of ISW, we plot in Fig. 4 (left) the ISW transform 
along the 60° latitude profile, called a spectrogram, for the satellite gravity disturbance 
(Nerem et al, 1994). We choose this small circle because it is approximately the latitude 
for both centers of Laurentide and Fenoscandia ice sheets. Thus, we can localize both 
major ice sheets simultaneously. As we can see from Fig. 4 (left) none of the small scale 
features in the gravity profile along the 60° latitude are lost in the spectrogram. The half 
angle of support in degree annotated on the left of the spectrogram characterizes the size 
or scale of the windows, the center wave number on the right indicates the central 
harmonic of the band-pass filter representing its bandwidth of a particular window. In 
other words, the spectrogram simultaneously reveals the spatial and harmonic 
localizations. Two major gravity lows are easily identified in Fig 4 (left): West Siberian 
Plain and Hudson Bay. The spectrogram in the West Siberian Plain area is rather 
dispersed; in contrast, the gravity low in the Hudson Bay area is heavily localized at 
harmonic degree 9, remarkably similar to the surface distribution of currently available 
models of the Laurentide ice sheet. This is one of the major arguments that led Simons & 
Hager (1997) to conclude that the gravity low in the Hudson Bay area does have a close 
association with incomplete rebound following the deglaciation. 

Since there is no direct observation of the history of the geoid off the shorelines, the 
only way we can reconstruct the RSL on the entire surface is through modeling. We 
adopt ice model ICE-1A as described above. 

Assuming the Earth is laterally homogeneous self-gravitating viscoelastic Maxwell 
body with its elastic parameters controlled by seismic inversion (Dziewonski & 
Anderson, 1981), the unknown radial viscosity profile is a major uncertainty if we want 
calculate the RSL at specified sites. But if our goal is to extract the long wavelength 
global pattern of geoid perturbation induced by deglaciation, the viscosity profile is much 
less a problem. In fact, we tried a number of viscosity models generated by our recently 
designed parameterization (Fang & Hager, 1999), and always come up with similar 
spectrogram as shown in Fig. 4 (right) where we deliberately leave the viscosity 
unspecified. The physics behind this observation is simple: the response of the solid 
Earth is laterally homogeneous if the rheology is laterally homogeneous. Thus, the 
heterogeneity of the geoid is mainly controlled by the irregular distribution of the ice 
sheets over the surface. Differences in viscous relaxation times resulting from 
differences in radial viscosity profiles only determines how fast detailed information of 
ice sheets is lost rather than altering the pattern revealed by the spectrogram. The ISW 
spectrogram of the predicted RSL in Fig. 4(right) is calculated at 9 kyr before present 
time (B.P.) A major portion of the ice sheets at this particular time has already melted. 
What is revealed in Fig. 4 (right) is primarily the viscoelastic relaxation of the solid 
Earth. We can clearly see the harmonic localization at about degree 9 in the central 
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region of the Laurentide ice, which is consistent with the previous result of Simons & 
Hager (1997), while the harmonic localization of the central region of the Fennoscandia 
ice sheet is at degree 16. The bandpass filters at the two specific centers are displayed in 
Fig. 5. The difference in harmonic localization is due to the fact that the Laurentide ice 
sheet covers a broader area than the Fennoscandian. Rebound signals at the center of an 
ice sheet reflect an average effect of the entire ice load. A broader spatial coverage is 
equivalent to a poorer spatial localization, which will result in a sharper harmonic 
localization, as controlled by the uncertainty principle. 

The significant bandwidth in the harmonic domain (Fig. 2) raises doubt whether the 
harmonic localizations revealed from the ISW spectrogram truly mean a single harmonic 



Fig. 4 (Left) observed satellite gravity disturbance along the 60 degree latitude and the wavelet 
spectrogram of the gravity disturbance. The translation operation of the wavelets are centered 
along the 60 degree latitude. Half angles of support (scaled windows) are annotated on the left, 
representing the spatial localization. The corresponding central harmonics are annotated on the 
right, representing the harmonic localization. The contour values are set for reference of 
sensitivity, they are not scaled to any physical unit. (Right) Predicted relative sea level (RSL) 9 
kyr before present time along the 60 degree latitude and the wavelet spectrogram of the RSL. 
Since the static sea level coincides with the geoid, we extend the meaning of sea level to the land 
as the geoid variation, thus, the IWS transform can be performed all over the globe. The 
calculation is based on the modified ice model, ICE-1A. We do not specify the viscosity here 
because for a laterally homogeneous model, the pattern of the wavelet spectrogram is insensitive 
to the viscosity profile. 
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dominant in the observed exponential like curves at the centers of the ice sheets. We 
look into this issue by comparing the decay time diagram of Mitrovica (1996) with the 
relaxation time diagram for single harmonics (Fig. 6). Mitrovica and Peltier (1993; 1995) 
parameterized the RSL curves in the central region of the major ice sheets by a single 
exponential form 

H(t) = A (exp(t/ T) - 1) (8) 


where H(t) denotes the height of the RSL at a specific site, A and T are the parameters 
determined by fitting the RSL curve at the site. The time, t, here is measured as positive 
backwards before the present time. Obviously, the RSL H(t) is a lumped signal of all 
harmonic degrees governed by the sea level equation (Farrell and Clark, 1976). To 
distinguish between the parameter T and the relaxation time for a single harmonic degree, 
x, we follow Mitrovica (1996), calling T a decay time. Assuming a two layer viscosity 
model with fixed lithosphere and ice model, one can find that the decay time is a function 
of the upper mantle viscosity r| u and the lower mantle viscosity, ri, 

T = T{r] um r} lm ) (9) 
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Fig. 5 Harmonic spectra of the IWS 
transforms at the centers of the two major 
ice sheets. Angle of supports (windows) are 
set about the same sizes of the two ice 
sheets. They all in the form of bandpass 
filter as discussed in the text. 
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Fig. 6. Relaxation diagram s of degree 16 
and 9. Shaded areas are adopted from 
Mitrovica (1996) for comparison. Thin dot 
lines are the relaxation diagrams for degree 
25 (top) and degree 4 (bottom). 
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For a fixed value of T, equation (2) determines an iso-contour in the rju—r), plane called 
the decay time diagram. In creating decay time diagrams, Mitrovica (1996) summed up 
harmonic degrees up to 256 to calculate the predicted RSL, H(t ), at the centers of disk ice 
loads simulating the Laurentide and Fennoscandia ice sheets. He then fit the observed 
RSL curves over the last 9 kyr at Angermanland Sweden, the center of the Fennoscandia 
ice sheet, and over the last 6.5 kyr at Richmond Gulf, Southeastern Hudson Bay, the 
center of the Laurentide ice sheet. 

For comparison, we generate a relaxation time diagram in Fig. 6 following a similar 
procedure to Mitrovica's (1996) decay time diagram: assuming a two layer viscosity 
model with a fixed lithosphere 

* = *0WL) (10) 

For a constant value of t, equation (10) determines an iso-contour in the plane, 
called a relaxation time diagram. 

We pick the degree 16 relaxation diagram to compare with the decay time diagram at 
Angerman land, and the degree 9 relaxation time diagram to compare with the decay time 
diagram at Richmond Gulf. Simultaneous matches in slope can be found in both the 
degree 16 (Angerman land) and degree 9 (Richmond Gulf) cases. Slopes are in fact the 
most diagnostic factor in identifying the harmonic degree in the relaxation diagram: the 
transition from lower degrees through intermediate degrees to higher degrees is very 
regular and can be clearly identified from the slopes of the iso-contours (compare the 
thin lines with the thick lines in Fig. 3). Thus, this simultaneous match in slopes strongly 
suggests that the decay times are indeed akin to the relaxation times of single harmonic 
degrees. 

In addition to the good fits in slope, the degree 16 relaxation diagrams (Fig. 6) match 
especially well with their Angermanland counterparts (Mitrovica's (1996) Fig. 4a) 
indicating a strong boost of degree 16 as a result of strong cancellations among the 
flanking degrees in the bandpass filter. As with the degree 9 relaxation diagrams 
compared to the Richmond Gulf decay diagrams, there are noticeable shifts in the 
relaxation diagram with reference to the decay diagram. For instance, the 6 kyr contour of 
the relaxation diagram locates approximately at the place of 6.6 kyr contour of the decay 
diagram, making about a 10% shift in iso-value. The largest of such shifts is about 20% 
by the 3 kyr contour. These discrepancies represent relatively larger ambiguity in relating 
the RSL to a single harmonic, as a result of poorer cancellation among the flanking 
degrees in the bandpass filter. 

Mitrovica and Peltier (1995) have found through numerical tests that if we only 
consider the time range in which the RSL change is dominated by free rebound, the decay 
time T is relatively insensitive to the details of the ice history. This observation can be 
easily explained if we interpret the decay times as strong signature of the relaxation times 
of single degrees, since the relaxation times are independent of ice models. Conversely, 
this observation indirectly supports our conclusion that the decay times, subject to 
ambiguities, represent single harmonic degrees. Our results indicate that the lateral 
heterogeneity in mantle rheology that affect the post glacial rebound is rather weak. 
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They also provide a new constraint on the relaxation spectrum first obtained by 
McConnell (1968) and rederived recently by Wieczerkowski et al. (1999) based on 
uplifted strandlines in the Fennoscandia area (Fig. 7). As shown in Fig. 7, the relaxation 
spectrum plays a crucial role in the inference of the mantle viscosity profile, especially in 
constraining the thickness of the lithosphere and the upper mantle viscosity structure. 
Detailed discussion on the relaxation spectrum and the viscosity structure is referred to 
Fang & Hager (2001). 


depth (km) 

0 1000 2000 3000 



Fig. 7. Relaxation spectra and the 
preferred viscosity models. McConnell’s 
(1968) spectrum with its error bars as 
interpreted by Cathles (1975) is plotted in 
the blue area. The yellow parallelogram 
represents a linear interpolation between 
degree 9 and degree 16 of our new results 
from wavelet analysis and the curve 
fitting based on the observed RSL at the 
centers of the ice sheets. The red line 
represents the spectrum from degree 10 to 
75 recently rederived by Wieczerkowski 
et al (1999). Plotted in the background as 
grey dots and inverse triangles are the 
viscosity model proposed by Simons and 
Hager (1997) and the relaxation spectrum 
associated with it. Viscosity r| is in the 
unit of 10**21 Pas. 


3.2 Results to be published 

The largest uncertainty in the GRACE data perhaps is the continental underground 
hydrological system. Its effect could span the entire normal lifetime (5 years) for the 
GRACE Mission. But the signatures of the hydrological system could be easily identified 
by a wavelet analysis, since the distribution of the underground reservoirs must be 
regionalized or localized within their residing continents. We would like to address the 
hydrology issue in a separate project. Currently are using a similar principle to quantify 
the ice mass balances in Greenland and Antarctic areas. We assume in the following 
discussions that secular terms have been removed from the time varying gravity field 
derived from GRACE. 

Our first objective is to identify the spatial and harmonic localization of the time 
varying gravity over Greenland and Antarctic. Note, spherical wavelets is a technique of 
2D to 4D mapping. The extra 2D space introduces redundancy for reconstruction. The 
redundancy could be overcome by what is called multi-resolution analysis (see below). 
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On the other hand, the extra 2D space provides room for zooming into small features and 
sharp edges like the ice sheet margin. A careful design of the scan contours could yield 
much more useful information. This idea of contour manipulation is already applied in 
choosing the 60° degree latitude small circle for characterizing the Laurentide and 
Fannoscandian Holocene ice sheets in Fig. 4 (right). For the Antarctic study, we choose 
another contour to characterize the two sub-domes of the Holocene ice sheet (Fig. 8). 
Obviously, these small circles (or great circles) are too simple for characterizing the 
migration and accumulation of the ice mass. 
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Fig. 8 Polar projections and wavelet spectrograms for predicted gravity disturbances due to last 
Quaternary ice sheets in the Antarctic area. The translation operation of the spectrogram is along 
the diagonal small circle indicated in the polar projections, in order to characterize the two small 
ice domes. Details of the settings of the spectrogram are referred to the caption of Fig. 4. 
Predictions are based on a Maxwell Earth with the viscosity model A in Fig. 1 and two ice sheets 
models. Ice-IA is specified in the figure 3 (see its caption). Ice-3G is from Tushingham & 
Peltier, 1992). Setting for the polar projections are the same as in Fig. 3. 
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We are now working on irregular contours for ISW scans. The most promising 
choices perhaps are the shoreline contours and iso-map contours. Since ice streams, 
including their surge and discharge, occur along the shorelines (e.g Hughes, 1998). The 
ISW spectrogram along the shoreline contours should be sensitive to the stream flow and 
discharge processes. This is important, since most of the uncertainties in the total 
imbalance come from the discharge processes (see the introduction). We are eager to see 
the results. But before we can apply the technique to analyze the GRACE data, a large 
amount of forward modeling has to be done in order to quantify the signature of stream 
flow and its resolution (the designed resolution for GRACE gravity is about harmonic 
degree 200, much less for its time variation therefore, resolution is a issue). We have to 
emphasize again that the ultimate resolution is limited by the raw data. Wavelets cannot 
create resolution for the raw data, it could only enhance and magnify the existing 
information blurred with conventional harmonic expansions. We use the digitized mass 
balance models (e. g. Reeh et al. 1999; Zwally & Giovinetto (2000); Vaughan et al. 1999) 
as the zero order input to simulate the GRACE gravity observations. We expect to derive 
the correlation coefficients between the mass change and gravity change in the wavelet 
transform domain. In order to investigate the sensitivity and resolution, we would 
manipulate the zero order mass balance model by artificially changing the discharged 
mass along the shorelines. An important part of this numerical experiment is with the ice 
shelves in the embayment confined by the shorelines and the calving fronts (the contour 
along which tabular icebergs are released into the open sea). We also conduct tests to see 
if the resolution for the volume of the stream flow could be improved by benchmarks 
with the local ground stations (see the local analysis). These works are time consuming, 
for the manipulations must be inspected by eye and digitized by hands. 

The iso-map contours can be generated either with the digitized topography or mass 
balance maps. Our ongoing works is similar to what just outlined for the shoreline 
contours. But here our focus is on quantifying gravity signatures of the accumulation and 
ablation processes for the ISW along the iso-map contours should be sensitive to 
accumulation and ablation. Melting is less volatile than surging and calving from ice 
stream flows, but its distribution is much more complicated. Most of the water bodies are 
under cover, large amount of meltwater are unable to reach the bed and flow out, because 
of regelation and refreezing (e.g. Lliboutry, 1998). Melting is more important a 
mechanism of discharge for Greenland than for Antarctic. Thus, we may regard the 
shoreline contours as the Antarctic contours, and the iso-map contours as the Greenland 
contours. While the ISW spectrograms from the two sets of contours are not directly 
comparable, the end results of inferred mass changes are comparable. Comparison 
between the two sets of results can serve as a test to see if our strategy really works. 

The second target of our on going effort is to quantify a specific mass balance 
signature identified from ISW transform by inverting it back to recover the original data. 
In wavelets terminology, it is an issue of data compression and reconstruction. The 
principle and general procedure of image compression is well known already (e.g. 
Resnikoff & Wells, 1998). For our purpose, we choose between two opposite quantizers: 
we either discard the specific mass signature to see its impact on the original data or to 
strip anything but the specific signatures to see how much the original data could recover. 
The recovery test is similar to the conventional harmonic analysis (e.g. James & Ivins, 
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1998; Bentely & Wahr, 1998) except in that the ISW signature in a specific location 
contributes significantly only to a selected band of harmonic degrees or bandpass filtered, 
while in a conventional harmonic analysis the mass signature contributes to the full 
harmonic spectrum. In other words, the ISW mass signature are characterized by its 
harmonic localization (see Fig. 8 for example). An alternative approach is to take some 
proxy quantities, such as the maximum of the local energy and the related local wave 
number (Bergeron et al, 1999; 2000). To complete the reconstruction, we have to go 
through the multresolution analysis (e.g., Freenden & Windheuser, 1996). The discovery 
of multiresolution analysis (Malat, 1989) perhaps is the ultimate reason for the popularity 
enjoyed by the wavelet technique today. The basic idea of multiresolution is to split the 
original signal f 0 into a “blurred” component/, with a coarser resolution and a “detailed” 
component d x with a regular resolution. The process repeats itself to split/, into f 2 and d 2 
. . .and go on and on. Similar recursive process can be formulated for the reconstruction. 
A remarkable thing is that the wavelets representing each d„ form a orthogonal basis with 
no redundancy. We skip the messy mathematical details except to say that, we have 
already completed the theoretical perpetration for a multi-resolution analysis based on 
our particular locally supported ISW, and are ready to write code for the recursive 
processes. The reason we would like to devote our energy to the multi-resolution 
analysis is that we could explore the ice mass balance problem from a different 
prospective, a global view. Our efforts in this part of the investigation are aimed at a 
multi-resolution inversion for mass balance (accumulation and ablation only) directly 
from simulated GRACE data. We mention the inversion problem just for a perspective; 
it is too premature at this stage to set a time frame for its completion. 

The price we pay with the ISW analysis is that the spatial localization for the mass 
balance signature is blurred, again, because of Heisenberg’s uncertainty principle. The 
mathematically exact spatial localization in a conventional harmonic analysis is just an 
illusion on the physical ground. It is exact only as we put it. But we do not know the 
exact distribution of subglacial water bodies, the drainage system, and so many other 
spatial distributions that associated with the mass balances. This is why we believe that 
wavelet analysis could make a difference in quantifying the total mass imbalances in 
Greenland and Antarctic areas: we trade the exactness in spatial localization, which we 
would suffer on the physical ground anyway, for a robust characterization of 
simultaneous spatial-harmonic localization. 

4. PROGRESSES OF GLOBAL ANALYSIS 

Time variation of lowermost harmonic degrees of the Earth’s gravity field, J nm , K nm can 
be regarded as benchmarks for the global mass redistribution (ice -i-ocean) in 
complementary to the local benchmarks of individual round measurements. In theory, 
GRACE is able to achieve better resolution than the satellite laser ranging (SLR). But the 
scheduled 5 year life time of GRACE already caused concerns that it may be too short for 
a complete separation between postglacial rebound and the mass migration of the 
underground hydrology system. A solution to this problem is to work together with SLR 
to extend its effective life (e.g. Wahr & Han, 1998). The SLR already has more than 25 
years of data collection and the effort is continuing. The low degree J nm , K n are sensitive 
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to mass redistributions not only on the surface, but also in the mantle and even in the core 
(e.g. Fang et al, 1996). The largest contribution from the interior is postglacial (e.g. Fig. 
9). But how sensitive the J nm , K n are to the model ice history is not fully understood, so 
we conducted, in collaboration with the CSR group at Texas U, the systematic 
calculations of the J nm , K n in the r| u ~r|, plane based on two ice models, ICE-1A and ICE- 
3G. We can see clearly in Fig. 9 that J„ is capable of distinguishing the two Holocene ice 
models. This is totally unexpected and encouraging. Recently, our collaborative efforts is 
directed to identifying the geophysical causes for J nm , K n . Our preliminary calculations, 
from published sources, of surface mass contributions to J n other than Greenland and 
Antarctic are listed in the following table. It is evident that Greenland & Antarctica are 
the major surface contributors. 



Satellite Obs. 

-2.7±0.4 

-1.2±0.5 

-1.1 ±1.0 

0.9±1 .0 

0.4±0.5 

-2.4±1 .4 

1.1±0. 

Mountain Glacier 

0.34 

0.25 

0.12 

0.25 

0.19 

0.24 

0.16 

Reservoirs 

-0.08 

-0.01 

0.11 

0.19 

0.20 

0.14 

<0.01 

Atmosphere 

0.17 

0.04 

0.18 

0.08 

0.08 

0.01 

0.01 

Ground water 

0.2 







Earthquake 

-0.022 

0.01 

-0.015 

-0.003 




Core rotation 

-0.04 


<0.001 


<0.001 




The SLR data processing and surface mass analysis are carried out by the Taxes 
group. Our job is to constrain the viscosity structure. Viscosity is a multidisciplinary 
issue, and unlikely to be resolved by the study of the J n . Nonetheless, it is important for 
us to understand how sensitive the J n are to the viscosity and how it may affect the 
estimates of the mass balance. Our preliminary results shown in Fig. 9 are based on a 
suit of two layer viscosity models. These two layer models are likely oversimplified, 
considering the real physical conditions in the Earth’s interior (e.g. Ranalli, 1998). We 
have come up with a new parameterization for the “realistic” viscosity models (Fang & 
Hager, 1999) which relates the models constructed from the microphysics of creeping to 
a set of parameters tangible with geophysical measurements. Two of such viscosity 
profiles are shown in Fig. 7. We are conducting similar calculations as shown in Fig. 9 
but with a suit of “realistic” viscosity models”. We expect to rule out some of the two 
layer viscosity models, if the predicted J n deviate significantly from those based on 
physical consideration and, at the same time, from those observed directly. We do not 
believe that such a procedure will enable us to sort out between the remaining two-layer 
viscosity models and the “realistic” viscosity models. Nonetheless, we could use the 
remaining two-layer mode to set the range of uncertainties for the estimates of the ice 
mass balance (the realistic viscosity profiles are a little difficult to handle for the purpose 
of setting error bars). 
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Ratio of lower-upper mantle viscosity 


Fig. 9 Time varying J n predicted base on postglacial rebound. Two ice models ICE1-A, and 
ICE3-G are used along with a suite of two-layer viscosity model with the lithosphere fixed at 120 
km and the upper mantle viscosity fixed at 10**21 Pas. The predicted are functions of the ratio of 
lower/upper mantle viscosity. Grey dots are from James & Ivins (1998). Shaded areas indicate 
the observed with their error bars from SLR (Cheng et al. 1997). 


Another issue we are investigating is how deep the J n can penetrate to sense the 
change of density associated with viscoelastic relaxation. It is more a geodynamic issue 
than geodetic. For the mass balance problem it allows us to explore the origin of the error 
source from the mantle. A method we plan to employ is the analysis of the relaxation 
time spectrum. As we can see from Fig. 7, the relaxation time spectrum is sensitive to the 
viscosity structure, and the lower degree relaxation times can penetrate deeper. We are 
trying to perturb different portions of the viscosity profiles from the elastic lithosphere, 
upper mantle down to 660km depth, and the lower mantle respectively to examine the 
responses of the relaxation time spectrum. We expect to provide through this experiment 
quantitative measures, in percentage, about the contributions from the lithosphere, upper 
mantle and lower mantles to the total error induced by mantle deformation in the 
estimates of the mass balances. 

Manuscripts for publication are prepared jointly and separately by the two groups. 
And the results shown in Fig. 9 have already been communicated with the geodynamic 
groups world wide. 

5. PROGRES OF SEA LEVEL AND TIDE GAUGE ANALYSIS 
5.1 Published results 

The sea level change associated with glacial isostatic adjustment is one of the classical 
problems in geophysics. As mentioned earlier, it is a coupled cryosphere-hydrosphere- 
atmosphere system. The melted ice mass spills into the sea, and rearranges the mass 
distribution on the surface of the solid Earth. The solid Earth in tum deforms in response 
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to the redistribution of the surface loads, and causes the sea level to adjust itself to reach 
an equilibrium with the remaining ice sheets and the deformed solid Earth. The melting 
of ice sheets has its astronomical causes and strong climatic consequences, while the 
deformation of the solid Earth depends crucially upon creep processes in the interior. 
Thus, the dynamics of the coupled system links a broad range of subjects, such as the 
evolution of planetary orbits, climatology, glaciology, mantle rheology, Earth rotation, 
and so on. 

The notion of a coupled system implies a global treatment, i.e. the volume of ice, sea 
water, and the solid Earth must be finite. Historically, the study of post glacial sea level 
started from a local approach (Haskell, 1935, 1937; Vening Meinesz, 1937; Niskannen, 
1948) for interpreting elevated beach traces in Fennoscandia. Locally, the ice load has 
forced the mantle rocks to flow laterally, causing considerable subsidence of the surface 
beneath the ice. In addition to this dynamic isostasy, the Earth and the ocean in a coupled 
system have to readjust themselves in responses to the redistribution of the ice load to 
maintain the toal gravitational potential energy of the system at a minimum. As it turns 
out, the restoring self-gravitation force against the Earth's readjustment plays an 
important role in the surface rebound (e.g. Peltier, 1989). 

The formulation for global analyses of post glacial sea level was begun by Cathles 
(1975), and completed by Farrell & Clark (1976) by the establishment of the so called sea 
level equation. Aside from the melting of ice and the responses of the finite solid Earth, 
the self-gravitation of the perturbed water load, i.e. the new sea level relative to the 
deformed sea floor, is also considered in the sea level equation to contribute to the 
formation of the new sea level itself. This self forcing mechanism makes the sea level 
equation non-linear. Approximations are implemented for solving the sea level equation 
at a particular time at a particular site (Farrell & Clark, 1976; Wu & Peltier, 1983; 
Nakiboglu et al, 1983). This situation is consistent with the fact that the observed sea 
level variationa are available only at a number of sites along continental shorelines and 
ocean islands. An effort has been made by Tushingham & Peltier, (1992) to compile the 
observed relative sea levels from various sources at more than 400 sites. Global solutions 
are proposed by Mitrovica & Peltier (1991) and Johnston (1993). 

A limitation of the site by site analysis is that it is difficult to reveal the physics of the 
coupled system by looking at individual sites. Often information about the viscosity and 
ice history is either too complicated to handle or too weak to use at individual sites. An 
example is the construction of the ICE-3G model (Tushingham & Peltier, 1991, 1992). 
ICE-3G is constructed by fitting the observed RSL curves with the predicted for a fixed 
viscosity structure at each individual site. Only 192 sites among nearly 400 were 
employed because the remaining 200 far field sites are not sensitive to the ice history. 
Among the 192 used, a large number of the sites located on the margin of former ice 
sheets are extremely sensitive to the uncertainties in both ice history and the viscosity 
(Nakada & Lambeck, 1987; Fang & Hager, 1996), hence the effects of viscosity and ice 
history at these sites are inseparable. It appears that substantial improvement can be 
made only by simultaneous constraints for both viscosity and ice history. A crucial step 
toward this goal is to obtain a global coverage of the geoid. 

Satellite altimetry, SLR, and especially the GRACE mission provide global coverage 
for time variations in geoid and sea level. But the satellite data are too recent to help 
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recreating the RSL history for the past 20,000 years. Isotope dating of the marine 
deposits along the elevated or subsided shorelines has been the major source for RSL 
data. Even in an ideal situation where the isotope data are collected continuously along 
the beaches all over the world, the data coverage is only one dimensional. Can we 
recover a 2-D distribution out of a set of I-D data? An issue relating to the answer is the 
smoothness of the 2-D field. The importance of smoothness to the recovery of the 2-D 
RSL is demonstrated by the following two extremes A perfectly smooth sea surface is the 
eustatic sea level. In this case the 2-D RSL can be recovered by the knowledge at only 
one site. A perfect non-smooth sea level is an extremely localized random type of 
distribution where the characteristic length for causal correlation between two points is 
nearly zero. The recovery of a total 2-D field from any kind of partial knowledge, 1-D, 
or 2-D is impossible. 

In this work we further explore the physics of the cryosphere-hydrosphere- 
atmosphere system. We first derived the sea level equation based on least potential 
energy principle. A global solution is also derived following the same spirit. Using the 
global solution, we next examine the smoothness of the global RSL. Our method is 
forward modeling on a laterally homogeneous Maxwell viscoelastic Earth using a 
number of admissible "realistic" viscosity models. We demonstrated that the horizontal 
smoothness of the RSL also relates to the sensitivity of RSL to the intra lithospheric 
structure. Parameterization of realistic viscosity models in itself an important issue. We 
propose a phenomenological parameterization. Based on microphysical considerations 
and seismic observations. This new parameterization can be directly constrained with 
geological observations. A lithosphere by definition suppresses viscous relaxation. We 
examine the effects of the lithosphere thickness on the sensitivity of RSL to "realistic" 
viscosity structures by comparing the predicted RSL for different lithosphere models. At 
a number of individual sites. We find that satisfactory convergence of global sea level 
solution can be reached at about harmonic degree 50. And the convergence appears 
independent of the ice model. This relatively lower tolerable truncation level is a 
consequence of global nature of the ice-sea-earth system. Variation of lithosphere 
thickness does not alter our previous conclusions (Fang & Hager, 1996) that there is a 
correlation of the RSL sensitivities between ice history and viscosity structures: at sites 
less sensitive to the ice model, the resolving power for viscosity structure is also less. 
Furthermore, models having a thicker lithosphere tend to permit better resolution of 
lower mantle viscosity. 

5.2 Results to be published 

The total ice mass imbalance from the continent is complementary to the total water 
imbalance in the ocean. The surface of an isothermal ocean at equilibrium can be derived 
based on the least potential energy principle (Fang & Hager, 1999). This steady state 
ocean is implicitly assumed in most of the geodetic sea level problems. The flowing 
ocean should also satisfy some least energy principle. The difficulty is how to specify 
those energies, especially the energies associated with the thermal transport processes. 
Dynamic effects superimposed on the steady ocean and local tectonics perhaps are the 
major error sources for the estimates of sea level variation from tidal gauge analysis. 
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We conducted a comparative study on the tidal gauge records from different regions. 
Our focus is on the Far-East coastal areas where the post glacial rebound effects are fairly 
small (see Fang & Hager, 1999 for specific sites in the Tokyo Bay area). A large number 
of tidal gauge records have been collected by the geodetic group in the Hong Kong 
Polytech. University over the last two decade (Ding et al., 2001). Fig. 10 shows the sea 
level variation derived from tide gauge records from the 1950’s to present at the Shanghai 
and Hong Kong stations. They both exhibit increases at about 2mm/yr. level, but the 
discrepancy, 0.6 mm/yr. is too conspicuous to ignore. This part of the investigation is in 
collaboration with the geodetic group at the Hong Kong Polytech. U. They have 
independent resources for their research projects under the International Program of 
“Asia Pacific Space Geodynamics”. 

The Hong Kong group has worked on the selection processes for the reliable sites and 
processing the data. Our expectation is to see if there are regularities from the 
discrepancies through the suit of sites. We are running our sea level code developed 
during this project to provide an improved steady state sea level. The present day steady 
state sea level have been calculated repeatedly (e.g. Conrad & Hager, 1997; Mitrovica et 
al, 2001) based on the equal potential theory of Farrell & Clark (1976) and simplified 
models of ice mass imbalances. We will use our least potential energy formulation (Fang 
& Hager, 1999) and new models of mass imbalance estimated based on cross constraints 
form local, regional, and global analysis. For a decadal time scale, an elastic reference 
Earth is adequate, thus, the computations are much simpler than for postglacial sea levels 
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Fig. 10 Upper panel, the 
monthly mean sea level 
change in Shanghai sea area 
during the time span of 

1945.0 ~ 2001.0 recorded 
from Wusong tide gauge 
station. The tide gauge data 
has been adjusted for land 
settlement by levelling 
measurements. The blue line 
shows a rising tendency with 
a rate of 1.79±0.19 mm/yr . 
Lower panel, Same as upper 
except that it is of Hong 
Kong sea area during the 
time span of 1954.0 ~ 

2001.0 recorded from North 
Point and Quarry Bay. The 
blue line shows a rising 
tendency with a rate of 

2. 3 9±0. 20 mm/yr. Source 
form Ding et al (2001) and 
Dawei Zheng of Hong Kong 
Polytech. U. 
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where an Maxwell viscoelastic Earth must be employed. We are waiting for the number 
to conduct a preliminary assessment of how reliable these tidal gauge results are in terms 
of relating them to the eustatic sea level change and mass imbalances. Results of this 
collaborative work will be published soon both jointly and separately by our two groups. 
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